5 de Mayo de 2022

Objetivo

Entender cómo usar los datos geolocalizados para problemas de marketing

En especial

  • Explicar la demografía de los visitantes de cada centro comercial

  • Usar los datos de movimientos de la gente para caracterizar la atracción de cada centro comercial

  • Estudiar problemas de logística/cierre centros comerciales

Qué es geomarketing

Recordamos: geomarketing es un tipo de “location intelligence” que, utilizando datos geolocalizados de tiendas, clientes y competencia permite:

  • Optimizar la inversión en acciones de marketing
  • Mayor conocimiento del mercado y de su segmentación geográfica
  • Optimización de la gestión espacial de recursos (tiendas, rutas, etc.)
  • Determinación de áreas de influencia de tiendas/empresas
  • Ayudar a las decisiones de inversión, compra o cierre de tiendas/oficias, etc.

Taller

El objetivo de este Taller es encontrar las zonas de atracción en los barrios de Madrid para cada uno de los Ikeas (u otros centros comerciales)

Para ello:

  • Vamos a utilizar datos geolocalizados de Twitter para determinar donde viven los usuarios que van a cada uno de los centros de Ikea de Madrid

  • Hacemos un modelo (fit) para determinar mediante el modelo de Huff cual es la probabilidad \(p_{ij}\) de que un usuario en un barrio \(i\) visite el centro de Ikea \(j\)

  • Con ese modelo determinamos las zonas de atracción de cada Ikea, asignando cada barrio \(i\) a un Ikea \(j\) si la probabilidad \(p_{ij}\) es mayor de las \(p_{i*}\)

Taller: Cargamos los datos geolocalizados

Vamos a utilizar una base de datos de tweets geolocalizados durante 2 años

load("./data/tweets_madrid.RData")
head(tweets_madrid)
##         lat       lon user                time
## 21 40.49162 -3.689221    1 2012-11-30 19:03:40
## 22 40.41296 -3.723460    2 2012-11-30 19:05:52
## 36 40.41916 -3.706140    3 2012-11-30 19:03:07
## 48 40.39049 -3.725962    4 2012-11-30 19:04:01
## 54 40.42712 -3.673200    5 2012-11-30 19:07:19
## 63 40.40596 -3.683690    6 2012-11-30 19:04:22

Taller: y el mapa de los barrios de Madrid

Cargamos los shapefiles de los barrios de Madrid

map_barrios <- read_sf("./data/Barrios Madrid/200001469.shp")
map_barrios$DESBDT <- iconv(map_barrios$DESBDT,from="Windows-1252", to="UTF-8")

Los transformamos a la proyección WSG84

mapb <- st_transform(map_barrios,4326)

Taller: y el mapa de los barrios de Madrid

Los visualizamos

leaflet(mapb) %>% addTiles() %>% addPolygons(weight=1)

Taller: poblacion de los barrios

También necesitaremos la población de los barrios de Madrid
http://www-2.munimadrid.es/TSE6/control/seleccionDatosBarrio

pob_barrios <- read.csv("./data/Madrid-Poblacion_Barrios.csv",
                      colClasses=c("character","character","integer",
                                   "character","character"))

Añadimos a la tabla del mapa una columna que sea la población (uniendo por el GEOCODIGO)

mapb <- merge(mapb,pob_barrios,by="GEOCODIGO")
head(mapb)
## Simple feature collection with 6 features and 7 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -3.722965 ymin: 40.40503 xmax: -3.690376 ymax: 40.43061
## Geodetic CRS:  WGS 84
##   GEOCODIGO CODBDT          DESBDT Distrito      Barrio Población     Fecha
## 1    079011 797664     011 Palacio   CENTRO     PALACIO     22471 10/1/2013
## 2    079012 797665 012 Embajadores   CENTRO EMBAJADORES     46264 10/1/2013
## 3    079013 797666      013 Cortes   CENTRO      CORTES     10717 10/1/2013
## 4    079014 797667    014 Justicia   CENTRO    JUSTICIA     16847 10/1/2013
## 5    079015 797668 015 Universidad   CENTRO UNIVERSIDAD     31714 10/1/2013
## 6    079016 797669         016 Sol   CENTRO         SOL      7577 10/1/2013
##                         geometry
## 1 POLYGON ((-3.709384 40.4224...
## 2 POLYGON ((-3.702836 40.4139...
## 3 POLYGON ((-3.696961 40.4189...
## 4 POLYGON ((-3.699416 40.4284...
## 5 POLYGON ((-3.712235 40.4302...
## 6 POLYGON ((-3.70416 40.42019...

Taller: y los datos de los centros comerciales

En concreto los centros comerciales en la comunidad de Madrid (http://www.madrid.org/nomecalles/)

centros <- read.table("./data/centros_comerciales.csv",header=T,stringsAsFactors = F)
head(centros,2)
##   CMUN     ETIQUETA         MUNICIPIO
## 1    5 Leroy Merlin Alcalá de Henares
## 2    7 Leroy Merlin          Alcorcón
##                                                                          DIRECCIÓN
## 1        Ctra Nacional II, km. 34 - Calle Tamajón s/n (Centro Comercial La Dehesa)
## 2 Avenida de Europa, 26 (Centro Comercial Parque Oeste - Ctra. Nacional V, km. 14)
##         lon      lat
## 1 -3.327188 40.50292
## 2 -3.849481 40.34330

Taller: determinación home usuarios

Estimamos donde está el domicilio de cada usuario a nivel de barrio: barrio más frecuentado para hacer los tweets.

Para ello añadimos una columna GEOCODIGO en la tabla de tweets de cada barrio:

points_tweets <- st_as_sf(tweets_madrid,coords=c("lon","lat"),crs=4326)
id_polygon <- st_within(points_tweets,mapb)
tweets_madrid$GEOCODIGO <- mapb$GEOCODIGO[as.numeric(id_polygon)]
head(tweets_madrid)

Taller: determinación home usuarios

Si va muy lento podemos utilizar también el paquete sp

WGS84 <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
points_tweets <- cbind(tweets_madrid$lon,tweets_madrid$lat)
pps1 <- SpatialPoints(points_tweets,proj4str=WGS84)
oo1 <- over(pps1,as(as_Spatial(mapb),"SpatialPolygons"))
tweets_madrid$GEOCODIGO <- mapb$GEOCODIGO[oo1]
head(tweets_madrid)

O podemos utilizar los tweets que ya están procesados en

load("./data/tweets_madrid_geo.RData")

Taller: determinación home usuarios

Nos quedamos solo con los tweets que suceden en alguno de los barrios

tws <- tweets_madrid[!is.na(tweets_madrid$GEOCODIGO),]

Es posible que haya usuarios que han venido a Madrid de vacaciones y hayan hecho muchos tweets. Nos quedamos solo con los tweets de aquellos usuarios que al menos haya estado en Madrid más de 5 meses. También nos quedamos sólo con aquellos que hayan hecho más de 10 tweets y menos de 10000.

tws$month <- month(tws$time)
users_table <- tws %>% 
                group_by(user) %>%
                summarize(ntot=length(user),nmonths=length(unique(month)))
users_sel <- users_table$user[users_table$nmonths > 5 & users_table$ntot > 10 & users_table$ntot < 10000]
tws <- tws[tws$user %in% users_sel,]

Taller: determinación home usuarios

Finalmente calculamos cual ha sido el barrio más frecuentado y cuánto ha sido frecuentado. Para cada usuario obtenemos cuantos tweets tiene en cada barrio

user_barrio_tweets <- tws %>% group_by(user,GEOCODIGO) %>% summarize(nt = length(user))

Ordenamos y calculamos el barrio más visitado

user_barrio_tweets <- user_barrio_tweets[order(user_barrio_tweets$user,user_barrio_tweets$nt,decreasing=T),]
user_barrio_most <- user_barrio_tweets %>% 
  group_by(user) %>% 
  summarize(ntot = sum(nt),
            nmax=head(nt,1),
            bmax=head(GEOCODIGO,1),
            nbarrios=length(unique(GEOCODIGO)))

Taller: determinación home usuarios

Solo consideramos aquellos en los que más del 30% de los tweets suceden en el barrio más visitado

user_barrio_most <- filter(user_barrio_most,ntot < 10000 & ntot > 10 & nmax/ntot > 0.3 & nbarrios > 3)

Taller: determinación home usuarios

Vamos a ver cómo se comparan el número de homes con la población de cada barrio

aa <- user_barrio_most %>% group_by(bmax) %>% summarise(nusers=length(user)) %>%
  ungroup()
merge(aa,pob_barrios,by.x="bmax",by.y="GEOCODIGO") %>%
  ggplot(aes(x=nusers,y=Población)) + geom_point() +
  xlab("Twitter users") + ylab("Población")

Taller: Estimación de áreas de atracción

Estimamos cuantos usuarios han ido a un centro comercial, desde que barrio y a qué distancia.

Vamos a probar con los de Ikea. Con el de Vallecas

centros_sel <- centros[centros$ETIQUETA=="Ikea",]
points_IKEA <- st_as_sf(centros_sel,coords=c("lon","lat"),crs=4326)

Vemos qué usuarios han geolocalizado un tweet a menos de 500 metros del centro de IKEA

points_tweets <- st_as_sf(tws,coords=c("lon","lat"),crs=4326)
distancia_a_Ikea <- st_distance(points_tweets,points_IKEA)
users_Ikea <- unique(tws$user[as.numeric(distancia_a_Ikea[,3]) < 500])

Taller: Estimación de áreas de atracción

Hacemos una tabla de los barrios de donde vienen esos usuarios y la distancia (al centro del barrio)

users_barrio_Ikea <- user_barrio_most[user_barrio_most$user %in% users_Ikea,]
agg_barrio_Ikea <- users_barrio_Ikea %>% group_by(bmax) %>% summarise(nn=length(user))

Juntamos estos datos con el mapa

mapb_Ikea <- mapb
mapb_Ikea <- merge(mapb_Ikea,agg_barrio_Ikea,by.x="GEOCODIGO",by.y="bmax")

Pintamos un choropleth con estos datos

mapp <- mapb_Ikea[!is.na(mapb_Ikea$nn),]
pal <- colorNumeric("magma",domain=mapp$nn)
m <- leaflet(mapp) %>% addTiles() %>%
  addPolygons(fillColor= ~ pal(nn),
              weight=1,fillOpacity = 0.8,
              popup=~paste(DESBDT,"<br> nusers = ",nn,sep=""))

Taller: Estimación de áreas de atracción

Veamos como es este choropleth

m

Ejercicio

Ejercicio 1:

  • ¿Desde qué barrio van más personas al centro de Ikea del norte?

  • ¿Cuál es la distancia media que recorren los clientes para ir a cada uno de los Ikeas?

Taller: ley de probabilidad de Huff

Finalmente, vamos a modelizar la atracción a cada uno de los centros comerciales (Ikea). Para ello vamos a utilizar la ley de Huff que dice que la probabilidad de que una persona que vive en \(i\) visite un centro comercial en \(j\) viene dada por

\[ p_{ij} \sim \frac{A_j}{d_{ij}^\beta} \]

donde \(A_j\) es lo atractivo que es el centro comercial \(j\) (normalmente se toma el tamaño), mientras que \(d_{ij}\) es la distancia del cliente \(i\) al centro \(j\).

Taller: ley de probabilidad de Huff

En nuestro caso vamos a empezar por ver cómo la distancia determina el número de usuarios que va a un centro comercial, mediante la ley

\[N_{ij} \sim \frac{1}{d_{ij}^\beta}\]

Para ello calculamos primero la distancia del Ikea al centro de cada barrio

cc <- st_centroid(mapp)
dd <- as.numeric(st_distance(cc,points_IKEA)[,3])

Taller: ley de probabilidad de Huff

Hacemos el fit (en logaritmos ya que es una función no lineal)

fit <- lm(log(nn) ~  log(dd), data=mapp)
summary(fit)
## 
## Call:
## lm(formula = log(nn) ~ log(dd), data = mapp)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.54144 -0.53346 -0.01751  0.50824  1.61527 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  15.4195     1.5635   9.862  < 2e-16 ***
## log(dd)      -1.5286     0.1705  -8.963 1.21e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6921 on 106 degrees of freedom
## Multiple R-squared:  0.4311, Adjusted R-squared:  0.4258 
## F-statistic: 80.33 on 1 and 106 DF,  p-value: 1.207e-14

donde vemos que \(\beta \approx\) 1.53

Taller: ley de probabilidad de Huff

Pero el numero de personas también depende de la población del barrio, asi que vamos a suponer una ley de la forma

\[N_{ij} \sim \frac{P_i^\alpha}{d_{ij}^\beta}\]

con lo que el modelo sería

fit1 <- lm(log(nn) ~  log(Población) + log(dd), data=mapp)
summary(fit1)
## 
## Call:
## lm(formula = log(nn) ~ log(Población) + log(dd), data = mapp)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.48464 -0.48147  0.01948  0.45129  1.65638 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    11.35004    1.77123   6.408 4.27e-09 ***
## log(Población)  0.39113    0.09633   4.060 9.46e-05 ***
## log(dd)        -1.51331    0.15935  -9.497 8.26e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6465 on 105 degrees of freedom
## Multiple R-squared:  0.5083, Adjusted R-squared:  0.499 
## F-statistic: 54.28 on 2 and 105 DF,  p-value: < 2.2e-16

Taller: ley de probabilidad de Huff

Dado que el fit es a un número de usuarios (variable discreta), probamos un GLM poissoniano con link logaritmico

\[\log (\lambda_{ij}) \sim \alpha log(P_i) - \beta log(d_{ij})\]

fit2 <- glm(nn ~ log(Población) + log(dd),data=mapp,family=poisson(link="log"))
summary(fit2)
## 
## Call:
## glm(formula = nn ~ log(Población) + log(dd), family = poisson(link = "log"), 
##     data = mapp)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -5.1388  -1.1508  -0.4280   0.5194   4.7271  
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)     8.85963    1.01115   8.762   <2e-16 ***
## log(Población)  0.76976    0.07101  10.841   <2e-16 ***
## log(dd)        -1.64339    0.06243 -26.324   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 1065.8  on 107  degrees of freedom
## Residual deviance:  246.5  on 105  degrees of freedom
## AIC: 610.31
## 
## Number of Fisher Scoring iterations: 5

Taller: Estimación de las zonas de atracción

Los resultados anteriores nos dan una forma funcional para el número de usuarios del barrio \(i\) que visitan el centro comercial \(j\):

\[N_{ij} \sim \frac{P_i^\alpha}{d_{ij}^\beta}\].

Vamos a utilizarla para determinar las zonas de atracción suponiendo que esa ley vale para todos los centros comerciales

  • Definimos una función para el modelo de gravedad de Huff
nij_Huff_model <- function(pts1,pts2,pob,alpha,beta){
  dd <- st_distance(pts1,pts2)
  prob_Huff_model <- pob^alpha/as.numeric(dd)^beta 
  prob_Huff_model
}

Taller: Estimación de las zonas de atracción

  • Calculamos \(N_{ij}\) para cada barrio y centro de Ikea
coord_barrios <- st_centroid(mapb)
nij <- data.frame()
for(i in 1:nrow(points_IKEA)){
  nij <- rbind(nij,nij_Huff_model(coord_barrios,points_IKEA[i,],
                                   mapb$Población,0.31,1.44))
}
colnames(nij) <- mapb$GEOCODIGO
head(nij[,1:8])
##         079011       079012       079013       079014       079015       079016
## 1 2.379124e-05 2.838968e-05 1.680126e-05 1.840339e-05 2.366358e-05 1.570227e-05
## 2 1.805954e-05 2.234352e-05 1.518181e-05 1.890398e-05 2.249597e-05 1.355679e-05
## 3 3.284778e-05 4.815319e-05 3.127393e-05 3.351510e-05 3.636477e-05 2.544260e-05
##         079021       079022
## 1 2.622335e-05 2.870612e-05
## 2 1.666106e-05 1.924022e-05
## 3 3.336490e-05 4.480561e-05

Taller: Estimación de las zonas de atracción

Para cada barrio, nos quedamos con el centro que da más número de visitantes. Eso define la zona de atracción. Y la visualizamos

centro_most <- apply(nij,2,which.max)
mapb$centro_most <- centro_most
factpal <- colorFactor(topo.colors(3), as.factor(mapb$centro_most))
leaflet(mapb) %>% addTiles() %>% addPolygons(weight=2,color = ~factpal(as.factor(centro_most)))

Ejercicio

  • Repetir el análisis de las zonas de atracción suponiendo que el IKEA de Alcorcón cierra.

  • Repetir el análisis para los centros de Leroy Merlin